int lbfgs(OptimFunc of, int n, double*xk, double*x0, int max ){
    int i, k, m, kk;
    double xkprev[n];
    double gk[n];
    double gkprev[n];
    double alfa = EPS;  // dlugosc kroku
    double *h[n];     // hesian
    double s[n];        // sk = xk+1 - xk
    double y[n];        // yk = grad(fk+1) - grad(fk)
    double *v[n];
    double *pkyksk[n];
    double *hl[n];
    double *h0k[n];
    double gammaK;
    int stop = 0;
    Vector Hkdeltafk;
    Vector pk;          // kierunek poszukiwan
    VectorHistory sk;
    VectorHistory yk;
    vhInit(&sk);
    vhInit(&yk);

    for (i=0;i<n;++i){
        h[i] = (double*)malloc(sizeof(double) * n);
        v[i] = (double*)malloc(sizeof(double) * n);
        pkyksk[i] = (double*)malloc(sizeof(double) * n);
        hl[i] = (double*)malloc(sizeof(double) * n);
        h0k[i] = (double*)malloc(sizeof(double) * n);
    }

    m = M;
    k = 0;
    kk = 0;
    mov(xk, x0, n);

    do{
	//k = 0;
        // Choose H0k (for example using 9.6)
        gammaK = 1.0;
        if (k>0){
            gammaK = countGammaK(sk.history[0], yk.history[0], n);
        }
        initMatrixIdentity(h0k, n, gammaK);
//        printMatrix(h0k, n);

        // Compute pk <- -Hk delta fk from Algorithm 9.1
        twoLoopRecursion91(of, n, xk, k, h0k,
                          &sk, &yk, /*out*/ Hkdeltafk);
        movNeg(pk, Hkdeltafk, n);

	//gradient(of, n, xk, gk);
	//movNeg(pk, gk);

//        printV("pk", pk, n);
        ///printf("pochodna w kierunku pk: %g\n", directedG(of,n,xk,pk,0));
        gradient(of, n, xk, gk);
        movNeg(gk,gk,n);
//        printV("-gk", gk, n);
        ///printf("pochodna w kierunku -gradient: %g\n", directedG(of,n,xk,gk,0));

        // kryterium stopu
        if (directedG(of,n,xk,pk,0) >= -EPS){
            /*printf("\n\nSTOP LBFGS %d!\n\n", kk);
            printV("xk", xk, n);
            printf("f(xk) = %g\n", of(xk, n));
            stop = 1;
            break;*/
            if ( (directedG(of,n,xk,gk,0) < -1*EPS) ){

                mov(pk, gk, n);
                k=0;
            } else {
                /// stop
                /// xk - optimal
                //printf("\n\nSTOP %d!\n\n", kk);
                //printV("xk", xk, n);
                //printf("f(xk) = %g\n", of(xk, n));
                stop = 1;
                break;
            }
        }

        // Compute xk+1 <- xk + alfak*pk, where alfak is chosen to satisfy the Wolfe conditions
        //alfa = findAlfaBisection(of, n, xk, pk, 5, 200);
        alfa = wolfe(of, n, xk, pk, 5, 200);
        mov(xkprev, xk, n);
        for (i=0; i<n; ++i) xk[i] = xk[i] + alfa*pk[i];

        // auto discard old history vectors

        // Compute and save sk <- xk+1 - xk, yk=delta fk+1 - delta fk
        subVectors(xk, xkprev, n, s);
        vhPush(&sk, s);
        gradient(of, n, xk, gk);  // gk = gradient xk
        gradient(of, n, xkprev, gkprev);  // gkprev = gradient xkprev
        subVectors(gk, gkprev, n, y);
        vhPush(&yk, y);

        // k <- k + 1
        ++k;
        ++kk;

        // debug info
//        printV("xk", xk, n);
        ///printf("f(xk) = %g\n", of(xk, n));
//        printV("gk", gk, n);
//        printf("sk history:\n");
        //vhPrint(&sk);
        //printf("yk history:\n");
        //vhPrint(&yk);


        //system("pause");
    } while(!stop);


    for (i=0;i<n;++i){
        free(h0k[i]);
        free(hl[i]);
        free(pkyksk[i]);
        free(v[i]);
        free(h[i]);
    }
    return 0;
}
